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ABSTRACT 

We re-analyze the precision radial velocity (RV) observations of HD 160691 (/j Ara) by the Anglo- Australian 
Planet Search Team. The star is supposed to host two Jovian companions (HD160691b, HD160691c) in long- 
period orbits (~ 630 days and ~ 2500 days, respectively) and a hot-Neptune (HD 16069 Id) in ~ 9 days orbit. 
We perform a global search for the best fits in the orbital parameters space with a hybrid code employing the 
genetic algorithm and simplex method. The stability of Keplerian fits is verified with the A^-body model of the 
RV signal that takes into account the dynamical constraints (so called GAMP method). Our analysis reveals 
a signature of the fourth, yet unconfirmed, Jupiter-like planet HD160691e in ~ 307 days orbit. Overall, the 
global architecture of four-planet configuration recalls the Solar system. All companions of /j Ara move in 
quasi-circular orbits. The orbits of two inner Jovian planets are close to the 2:1 mean motion resonance. The 
alternative three-planet system involves two Jovian planets in eccentric orbits (e ~ 0.3), close to the 4:1 MMR, 
but it yields a significantly worse fit to the data. We also verify a hypothesis of the 1 : 1 MMR in the subsystem 
of two inner Jovian planets in the four-planet model. 

Subject headings: celestial mechanics, stellar dynamics — methods: numerical, A^-body simulations — plane- 
tary systems — stars: individual (HD160691) 



1. INTRODUCTION 

The star HD 160691 (/j Ara) is a Sun-like main-sequence 
dwarf monitored by the long-term, precision radial velocity 
(RV) surveys. It has been observed over more than 7 years by 
the Anglo-Australian Telescope (AAT) Planet Search Team 
and by the Geneva Planet Search Team (CORALIE and 
HARPS spectrometers). The work of the AAT team lead 
to the discovery of a J upiter like c ompanion HD160691b in 
about 630 days o rbit (iButler et al.l 12.0011) . One year later, 
iJones et all (120021) confirmed the Jovian planet and discov- 
ered a linear trend in the RV data revealing a signatur e of the 
second , more distant body. In the next paper, Mc Carthy et alJ 
(2004) published a new orbital solution with the orbital pe- 
riod of the long -period planet HD 160691c about 3000 days 
and th e eccentricity e c ~ 0.57. The same year, ISantos et al.l 
( 2004), using observations done with the ultra precise HARPS 
spectrometer, announced ~ 14 Earth-mass planet HD160691d 
in ~ 9 days orbit. That discovery is a breakthrough in the 
field as the long-term precision of spectrometers approaches 
1 m/s. Actually, the instrumental errors are much smaller 
than the RV variability (stella r jitter) induced b y the Sun- 
like stars themselves. Recently, Bu tler et alJ 12006) published 
a new updated set of 108 observations of /j Ara, spanning 
255 1 days (about 7.5 yr). Thanks to the updates in the UCLES 
instrument installe d at the AAT and the software pipeline 
jButler et al J 120061) . the long-term precision of the measure- 
ments is amazing. It also reaches 1 m/s at the end of the ob- 
servational window and is kept at the mean level of ~ 2.8 m/s 
over its whole length. 

This paper is devoted to the analysis of the RV observations 
of fi Ara using the so called GAMP approach (an acronym 
of genetic algorithm with MEGNO penalty) that makes ex- 
plicit the use of Newtonian, self-consistent A^-body model of 
the stars' reflex motion, dynamical properties of the plane- 
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tary system and the Copernican Principle ( Gozdz iewski et all 
2003, 2005). This algorithm makes it possible to derive mean- 
ingful bounds on the elements of the outermost planet in spite 
of the fact that the data cover only a part of the longest orbital 
period. Without the dynamical constrains, the kinematical as 
well as pure Newtonian best fits to the /j Ara RV tend to show 
large eccentricity of the outermost companion and then the 
system becomes catastrophically unstable. 

We verify the results of the previous papers based on a 
much smaller number of relatively less accurate RV obser- 
vations. The results of the analysis of the n ew RV data 
set gre atly improved and extended over time by Butl er et all 
(2006) lead us to a conclusion that the /j Ara may host four 
planets in quasi-circular orbits, including a new, Jupiter-like 
object in ~ 307 days orbit (HD160691e) 3 . The new best fit 
solution describes the orbital architecture of this extrasolar 
system as very different from the previous ones. 

Shortly after submittin g the manuscript , we learned about 
an independent work bv lPepe et al.l (1200 6) who announced a 
very similar orbital solution and also the fourth planet in the 

Ara system. These authors study a different set of the RV 
measurements but the conclusions of the two papers are in an 
excellent agreement. Still, in this paper we focus our atten- 
tion on the analysis of the AAT data along the lines of the 
originally submitted manuscript. However, the last Section 5 
is devoted to a preliminary study of the full availab le data 
set, inc luding the new measurements published bv lPepe et alJ 
(2006). In particular, we verify a hypothesis of the 1:1 MMR 
in the subsystem of two Jovian planets. 

2. FITTING MULTI-PLANET CONFIGURATIONS TO RV DATA 

According to the previous papers devot ed to u Ara, the time 
range of the updated data set published by B utler et ahl (l2006|) 
should be already close to the orbital period of the outermost 
companion. In that case, we can try to recover a good approx- 
imation of the system parameters by modeling the RV signal 

3 In this paper we use a widely adopted (but unofficial) convention for 
naming the extrasolar planets with lower-case Roman letters starting from 
"b", in the order of their announcement. 
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with the Keplerian orbits. Although for mutually interacting 
systems the kinematic model often leads to unstable orbital 
configurations, thanks to the numerical simplicity it can be 
very helpful to rapidly determine putative solutions and inter- 
esting ranges of orbital parameters. 

The contribution of every planet to the reflex motion of the 
parent star at time t is the following Smart ( 1949): 

V r (t) = #[cos(co + v(r))+ecosco]+V , (1) 

where K is the semi-amplitude, co the argument of pericenter, 
v(f) the true anomaly involving implicit dependence on the 
orbital period P and the time of periastron passage T p , e the 
eccentricity, and Vb is the velocity offset. Some argue that 
it is best to interpret the derived fit parameters (K,P, e,co, T p ) 
in terms of Keplerian eleme nts and minimal masses r elated 
to Jac obi coordinates tLee & P eale 2003; Gozdziews ki et alJ 
2003). We follow their reasoning when calculating the orbital 
elements from the primary fit parameters. 

The extensive exploration of the multi-parameter 
(Xy) 1 / 2 space is efficient enough if one applies a kind 
of a hybrid optimization (IGozdziewski & Konackil l2004t 
Gozdziewski & Migaszewski 2006). The idea of this algo- 
rithm relies on two steps, global and local ones. In the first 
step, we search for potentially good (but not very accurate) 
solutions in a global manner. During the second step these 
solutions become initial conditions to a precise and fast 
local algorithm. The single program run starts the genetic 
al gorithm (GAs) . In particular, we apply the PIKAIA code 
by Charbonneau ( 1995). The GAs have imp ortant advantage s 
over more popular gradient-type methods (Press et al. 1992). 
The power of GAs lies in their basically global nature, 
the requirement of knowing only the (x 2 ) 1//2 function, and 
the ease of constrained optimization; GAs permit defining 
parameter bounds according to specific requirements, or 
adding a penalty term to (% 2 ) 1//2 (Gozdziewski et al. 2006). 
However, the best fits found with GAs are (in principle) not 
very accurate in terms of (%l) 1 ^ 2 or rms, so we refine them 
with another non-gra dient algorithm, the simplex method 
of Melder and Nead (lPressetalJ lT992). Usually, we run 
such hybrid procedure thousands of times, and then we 
analyse the ensemble of gathered fits. That helps us to detect 
local minima of (x 2 ) 1//2 mat ^ sometimes very distant in 
the parameter space and to get reliable approximation of 
the glo bal topology of (y 2 ) 1 / 2 . We tested the code exten- 
sively (Gozdziewski & Migaszewski 2006), and we found 
many examples confirming its robustness and reliability. 
Remarkably, the code works with minimal requirements for 
user-supplied information: basically, one should only define 
the model function [so called fitness function, usually equal 
to 1/ (Xv) 1 ^ 2 ] — conveniently, it is the same for the GAs and 
simplex, and to determine (even very roughly) parameter 
bounds for the assumed number of planets. We underline that 
the code is FFT-free, and by its construction, it works without 
any a-priori determination of the orbital periods. 

In some cases (like strongly resonant or interacting sys- 
tems, noisy data, a small number of measurements) the 
kinematic fits may lead to unrealistic, rapidly disrupting 
configurations. Then a more elab orate TV-body Newtonian 
model of the RV should be a pplied ( Rivera & Lissauerll2001t 
Laugh lin & Chambersl200ll) . Yet the hybrid optimization can 
be still used as a general approach of exploring the parameter 
space (only the model function is changed). Actually, even 
more general modeling of the RV data relies on the elimi- 



nation of the unstable (strongly chaotic) solutions during the 
fit process. That self-consistent approach follows the Coper- 
nican Principle and the complex structure of the phas e space 
predic ted by the Kolmogorov-Arnold-Moser theorem (lArnoldl 
1978). The idea of the dynamical fits with stability constraints 
relies on modifying the (y 2 ) l ' 2 function by a penalty term em- 
ployin g an efficient fast indicator MEGNO (ICincotta & Simol 
2000). We describe that method (GAMP) in detail in our p ast 
works (Gozdziewski et al. 2005; Gozdziewski et al. 2006). 

To obtain reliable estimates of the fitted parameter er- 
rors, the internal measurement errors should be rescaled 
( iButler et alJ 12004 according to a 2 = + where a m 
and Oj is the internal error and adopted dispersion of stellar 
jitter, respectively, and a is the joint uncertainty. Typically, 
we choose a\ following the estimates for Sun-like dwarfs by 
Wright (2005), or we use the value adopted b y the discove ry 
teams. The jitter estimate for /j Ara by Butler et al. (2006) is 
3.5 m/s and we use this value in our calculations. The anal- 
ysis of the short time-scale RV variability of the parent star 
may be found in Bouchv et al] (120051). The physical proper- 
ties of the parent star ar e discussed in McCarthy et alJ (120041) 
and lSantos etail (120041) . 

3. THREE-PLANET MODEL OF THE RV 

In (IGozdziewski et al.1 l200l 120051) we carried out an 
exten sive analysis of the RV measur ements published by 
Uones et al] (120021) and lMcCarthv etalJ ( 12004 assuming that 
there exist two Jupiter companions of p Ara. By employing 
the dynamical TV-body model and GAMP, we obtained mean- 
ingful limits on the barely constrained orbital parameters of 
the putative outermost planet. According to our results, a c 
should be roughly greater that 4 AU and e c < 0.4 in the range 
of the smallest permissible s emi-major axes. Th e new pre- 
cision RV data published by Butl er et aTl (|2*006) give us an 
excellent opportunity to verify these conclusions. 

Figure \l\ shows the results of the hybrid search for the 
putative three-planet configurations, in terms of multi-planet 
Keplerian model, assuming that the innermost planet has 
Pd G [7,12] days and £ [0,0.3], and that two Jupiter- 
like planets are in orbits with Pb G [100, 1200] days, P c G 
[1200,8500] days, and e b , c G [0.0,0.8], respectively. These 
safe assumptions are consistent with the results of a few pa- 
pers devot ed to u Ara. In part icul ar, we rely on the ca reful 
analysis bv lSantos et al.1 ll2004 and lBouchv et al.1 J2005I) that 
revealed the hot-Neptune HD 16069 Id. We might expect that 
its signal K& ~ 4 m/s is comparable to the error level of the 
AAT data, nevertheless we decide to add this planet to the 
model, not to avoid the a-priori information on the system 
architecture. It appears that the hot-Neptune signal improves 
the rms by ~ 0.5 m/s, so it could be important to obtain a 
precise solution for the whole system. 

The best fits obtained in the search are illustrated by pro- 
jections onto the planes of particular parameters of the RV 
model, Eq. 1. Here, we choose the (P,K)- and (P,e) -planes. 
Marking the elements within the formal la, 2a, 3a confidence 
intervals of the best-fit solution (signed by two crossing lines), 
we have a convenient way of visualizing the shape of the local 
minima of (y 2 ) 1//2 an d obtaini ng realistic and reasonable esti - 
mates of the parameter errors (Bevington & Robinson 2003). 
Figure 1 illustrates the parameters of ~ 1000 different fits 
within the 3o confidence interval of the best Fit I (its param- 
eters are given in Table 1). Remarkably, the orbital period 
Pd ~ 9.637 days and the semi-amplitud e A"j ~ 3 m/s a re very 
close to the independent estimates by Santos et al. (2004), 
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on the basis of HARPS measurements. Thus in spite of the 
RV contribution of the hot-Neptune planet d (of the inferred 
~ 1 1 Earth-masses) being on the noise level ~ 4 m/s, the 
planet is already "visible" also in the AAT measurements. The 
two-planet Keplerian model of the RV yields 4 (x 2 ) '^ 2 ~ 1-1 1 
and an rms ~ 4.45 m/s, so the signal of HD160691d improves 
the fit by ~ 0.5 m/s. Yet in the next section, we are trying to 
find much better arguments supporting this claim. 

In turn, the elements of the Jupiter-like companions seem 
constrained very well. For a reference, the synthetic curve 
and the data points are illustrated in Fig. [3 An rms of Fit I is 
~ 4 m/s and its (x 2 ) 1 ^ 2 ~ 1- The best-fit three-pl anet solution 
found here is very similar to the one quoted by Butl er et ail 
(2006), see their Table 3. The orbital period ratio close to 
4: 1 suggests a proximity of the Jovian planets to the 4: 1 mean 
motion resonance (MMR). Unfortunately, the system is again 
catastrophically unstable due to a large eccentricity of the 
outer planet (e c ~ 0.47) and a proximity of both orbits to the 
collision zone, i.e., the area close to the planetary collision 
line determined by a\,{\ +<?b) = a c (l — e c ). 

However, it is well known that orbits involved in 
low-order MMRs m ay be stable even if they are 
crossing each other (Jietal. 2003; Beauge et al. 2003; 
Psvchovos & Hadiidemetriou 2005). In Gozdziewski et al. 
(2006) we re-analyse the RV of HD 108874 dVo gt et al.l2005l) 
that appears to host two Jovian planets very close to the same 
type of the 4:1 MMR. The dynamical map of this system 
reveals that the eccentricity of the outer planet could be as 
large as 0.7 in the stable resonance island, although already 
for e c ~ 0.4 the osculating orbits would cross. In the same 
way, the exact 4:1 MMR could explain large e c ~ 0.5 in the 
three-planet /u Ara system. To examine more carefully the 
stability of the best-fit configuration, we computed dynamical 
maps (Fig. [3} i n the (a c ,e c )-plane, in terms of the Spectral 
Number (SN) (Michtchenko & Ferraz-Meildl20"o"ll) and the 
maxe indicator (the maximal eccentricity attained during the 
integration time-span). In particular, every point in these 
maps represents an initial condition that was integrated over 
~ 10 5 yr (~ 10 4 P C ). The dynamical maps reveal a few 
dominant low-order MMRs: like 4b: lc, 9b:2c and 5b: lc 
in the neighborhood of Fit I, marked by a crossed circle 
(«bb : « c c means the «b : n c MMR of planets "b" and "c"). 
Clearly, the best-fit Keplerian configuration lies in a strongly 
chaotic zone 5 . 

A simple change of the Keplerian best-fit <? c to ~ 0.25, 
providing a stable system, leads to a significant increase 

of (Xv) 1,/2 ' so we tr i e d to an optimal stable solution 
with GAMP. Due to a large CPU requirement caused by 
the short-period orbit of the innermost low-mass HD 16069 Id 
that planet has been skipped in this test. We searched 
for the two-planet solutions only. In the penalty function 
( IGozdziewski et all 12003) we integrated the MEGNO over 
~ 1Q 3 P C . It is a relatively short time but it enables us to rule 
out strongly chaotic (and rapidly disrupting) systems. The 
results of that search are illustrated in Fig. @] were projec- 

4 The parameters {K,P.e,a>, T — 7b) in Eq. 1 are (37.52 m/s, 630.31 days, 
0.269, 259.63 deg, 13401.53 days), and (18.10 m/s, 2499.16 days, 0.466, 
184.04 deg, 11032.98 days), V = -0.03 m/s, T =JD2,440,000. 

5 In I Gozdziewski et al. 2005) we already derived a very similar solu- 
tion to the ones quoted by Butler et al. 1 2006 ) and found in this work, with 
a c ~ 3 .8 AU and large e c ~ 0.6, on the basis of RV data from McCarfhyjjLay 
1 2004) extended by measurements published graphically in San tos et alj 
1 2004). However, we could not find any stable solution in its close neigh- 
borhood. 



tions of the best-fit parameters onto the dynamical maps in 
the (a c ,e c )-plane are shown. Only the solutions within the la 
confidence interval of the best fit (marked by the largest cir- 
cle, see also caption to Fig.|4]for its osculating elements) are 
shown. Let us note that this best-fit solution has been refined 
by GAMP integrations over ~ 25 , 000P C and the solution ap- 
pears to be rigorously stable. An rms of these two-planet fits is 
~ 4.7 m/s and their (x 2 ) 1//2 ~ 1.17. The scatter of the best-fit 
parameters is small but we cannot decide whether the system 
is locked in the exact 4: 1 MMR. More likely it evolves close 
to its separatrix (outside the resonance island). Note that the 
fits with the largest e c (to the left of the resonance island) are 
in fact mildly chaotic. The presence of the Neptune-like com- 
panion does not lead to any qualitative changes in the dynam- 
ical character of the best fit configurations although it yields a 
lower rms ~ 4.4 m/s (for stable configurations). 

The results of modeling the RV data by the three-planet 
configurations seem in an o verall agreement with the conclu- 
sions of our previous work (Gozdziewski et al. 2005). How- 
ever, we found that the data published in McCarthy et al. 

(2004) rather exclude the possibility of a stable 4:1 MMR be- 
tween the Jovian planets. The acceptable (stable) solutions 
should have a c roughly not smaller than 4 AU. The corre- 
sponding orbital period is significantly longer, by ~ 500 days, 
from the current apparently very precise estimate of P c ~ 
2500 days found on the basis of the new data set. Still, al- 
though the observational window already covers about one 
outermost period P c and the data strongly constrain (x 2 ) 1//2 i 
both the kinematic as well as the Newtonian model of the RV 
yield catastrophically unstable orbital configurations. 

That lead us to look for an explanation of the strange in- 
consistency. The most natural one could follow from the exis- 
tence of an additional planet that has been hidden up till now 
due to the small number of measurements and their signifi- 
cant errors (~ 4 m/s, as quoted in the older papers by the AAT 
Tea m). The problem reminds us of the study of the HD 37124 
data dVogt et all2005l IGo zdziewski et al. 2006). For this star, 
the two-planet fits are strongly unstable due to the e xtreme ec- 
centric ity of the outer companion, ~ 0.7. Recentlv.lVo gt et all 

(2005) has shown that the assumption of three-planets makes 
it possible to improve the fits and, simultaneously, the best-fit 
orbits become close to circular one s. The system can also be 
easily stabilized in such a regime (IGozdziewski et alJ 2006) 
without any degrada tion (y 2 ) 1 / 2 and the rms. Further, we 
follow the results of iBouchv et al.l (120051) who detected the 
short-period planet d around /j Ara with ultra-precise HARPS 
observations. In particular, we are attracted by the analysis of 
the short-time scatter of RV during several subsequent nights 
(see their Table 1 and Fig. 1,2). That statistics measures the 
RV variability imposed by the star activity itself. The stan- 
dard deviation of these variations over one night is between 
1.5-2.5 m/s. It could mean that aj is in fact less than 3.5 m/s 
that we adopted here. Instead, assuming aj ~ 2 m/s we got 

(Xy) 1 / 2 ~ 1.4 for the best-fit three-planet model and the rms 
~ 4 m/s has ~ 1.2 m/s excess over (a) ~ 2.8 m/s (the mean 
of O m is ~ 1.9 m/s). In such a case, the three-planet model is 
not fully adequate for explaining the RV variability and that 
the presence of an additional, yet undetected planet, is statis- 
tically justified. 

4. FOUR-PLANET SYSTEM AND ITS ORBITAL STABILITY 

To verify the hypothesis about the four-planet system, we 
again used the hybrid code driven by the kinematic model of 
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the RV. This time we assumed that all orbital periods are in the 
range of ~ [8,6500] days. The statistics of the best fits gath- 
ered in the search are illustrated in Fig. [5] We marked ~ 1000 
different solutions yielding (x 2 ) 1//2 within the 3o confidence 
interval of the best fit with P c ~ 4000 days, (x 2 ) 1/2 = 0.626, 
and an rms of 2.276 m/s (marked by crossi ng lines). How- 
ever, after a check by MEGNO integrations ( Cincotta & Simo 
2000), we found that this fit leads to a chaotic configura- 
tion, so we selected another stable (quasi-periodic) solution 
with very similar (Xv) 1 ^ 2 = 0.627, an rms of 2.276 m/s and 
P c ~ 4500 days. Its orbital parameters are given in Table 2 and 
we call it Fit II from hereafter. We use that initial condition 
to demonstrate some orbital properties of the /j Ara system. 
We tried to refine this solution by the self-consistent A^-body 
model of the RV but we could not improve it significantly; the 
osculating elements also did not change. 

In general, the four-planet solutions have much smaller 
rms' (by ~ 1 .7 m/s) than the three-body configurations and 
they impose the new, ~ 0.5 mj-mass Jovian planet, much 
closer to the star than the companion HD160691b in ~ 
650 days orbit detected a few years ago. The most striking 
feature of the four-planet configurations is low eccentricity 
of all orbits. They are roughly less than 0.2 for all Jovian 
planets. The innermost hot-Neptune planet d has the initial 
ed ~ 0.2, nevertheless, it is not well constrained and any value 
in the (assumed) range [0,0.3] is equally likely in terms of the 
la confidence interval. Simultaneously, the orbital periods 
of the two inner Jupiter-like companions, Pb ~ 646 days and 
P e ~ 307 d appear to be bounded very well with the accu- 
racy range of a few days. This is not the case for the outer- 
most planet c — the la error of P c is about 1500 days. Thus 
the four-planet model changes completely the topology of 
(Xv) 1 ^ 2 - Yet the outermost planet would have the semi-major 
axis very s imilar to tha t of Ju piter. According to the con- 
clusions of Santos et al. (2004), the orbit of planet d should 
be almost edge-on. Hence assuming that the entire system is 
coplanar, we may expect that the minimal masses determined 
from the Keplerian fits are likely close to the real ones. 

Figure[6]illustrates in subsequent panels the synthetic curve 
of the best-fit system (Fit II, Table 2) and the period-phased 
RV signals of the planets d, e, b, and c. The resulting curve 
closely follows all the measurement points. The rms is only 
~ 2.3 m/s. As we expected, its value is comparable with 
the joint error a if we assume that Oj ~ 2 m/s. This indi- 
cates a statistically perfect solution. The next panel shows 
the raw Lomb-Scargle periodogram (Pr ess et aljJ l992) of the 
data. Besides the dominant signal of planet b, there are only 
visible peaks about 32 days and 225 days. Simultaneously, 
the periods of ~ 9 days and ~ 307 days seem to be com- 
pletely absent in the periodogram. The period of ~ 32 days 
is close to the rotational period P rot of the star deriv ed from 
the index log/^ K = -5.034 by iSantos etaP (120041) . How- 
ever, it is not consist ent with another estimate of ~ 22 days by 
iBouchv et al.l (|2005 ) who analyzed the short-term variability 
of the RV spectrum of /j Ara. We have no good explanation 
for the ~ 225 days period as we did not find any reasonable 
solution consistent with its value. Most likely, it is an alias 
of the ~ 32 days period. Instead, there exist other apparently 
precise fits having periods uncorrelated with the periodogram, 
for instance, (2.332 m/s, 90.19 days, 0.30, 3.739 deg, 
14391.686 days), (12.298 m/s, 517.578 days, 0.581, 
157.537 deg, 13841.816 days), (39.010 m/s, 621.323 days, 
0.267, 259.890 deg, 10908.948 days), and (21.580 m/s, 



2459.212 days, 0.475, 183.853 days, 13542.752 days), in 
terms of parameter tuples of Eq. 1, yielding an rms ~ 2.7 m/s 
(T p are shifted by JD2,440,000). However, these solutions are 
very unstable. We analyse such alternative solutions in Sec- 
tion 5.1. 

Another unusual coincidence of periods is the ratio P e /Pi ~ 
32. That may also indicate an aliasing in the data. Never- 
theless, the signal of K e ~ 14 m/s exceeds by a few times 
the errors of data and our experiments show that it cannot 
be fitted well by the signal of the hot-Neptune planet itself. 
Still, it remains possible that the periodic signal of the putative 
planet e is an alias of the rotational period because P e ~ 10P rot 
if Pot ~ 32 days and P e ~ 14P rot if P rot ~ 22 days. Then the 
RV variability could be attributed to a mo ving spot on the star 
surface. However, in the observations of Santos et al. (12001 
there is no sign of abnormally large variations of the RV (ex- 
ceeding the amplitude of the hot-Neptune signal) during about 
80 days, covering already a few P rot . It is a good argument 
justifying the planetary origin of the signal with the period of 
~ 307 days. 

Finally, we investigate whether the hot-Neptune 
HD160691d can be reliably detected in the AAT data 
alone. Because we fit a multi-period orbital model to a small 
number of measurements, there is always a risk of generating 
a periodic signal through random fluctuations. To determine 
the confidence level of a weak signal i n the data, we app ly the 
so called test of scrambled residuals (Butler et al. 2004), see 
also our paper Gozdziewski & Miaaszewski (2006). Having 
the best fit (Fit II) in hand, we remove the synthetic RV 
contributions of the Jovian planets from the measurements. 
Next, we randomly scramble the residuals, keeping the exact 
moments of observations, and we search for the best-fit 
elements of the putative planetary signal. To speed up the 
search, we limit the period range to [2,128] days. We use 
again the hybrid code. Thanks to a large population size 
(1024) and the number of generations (256), the GAs reliably 
find the best fits to the single-planet model; yet the fits are 
then refined with the simplex. After many repetitions of such 
a procedure, we get a Gaussian-like histogram of (x 2 ) 1 ^ 2 of 
the best fits to the synthetic sets of residuals. If the real data 
were uncorrelated, the value of their (x 2 ) 1,/2 should be found 
in the range spanned by the histogram. That histogram of 
~ 36,000 Keplerian fits to scrambled residuals is shown in 
the right-bottom panel of Fig. |6] Clearly, the likelihood that 
the residual signal is only a white noise is negligible. For a 
reference, we computed the Lomb-Scargle periodogram of 
the original residuals (the left-bottom panel in Fig. [6}. The 
strongest peak at ~ 9.637 days is consistent with the estimate 
of Pa in the hybrid search (see Fig.|5}. These results support 
independe ntly the de tection of the hot-Neptune HD 16069 Id 
by Santos et al. (2004) and prove the excellent quality of the 
AAT data. 

4.1. Stability analysis 

To investigate the long-term stability of the best-fit configu- 
ration, we employed the MEGNO indicator computed by the 
symplectic algorithm (Gozdziewski et al. 2005). In the first 
test, we integrated the entire four-planet system (with the ini- 
tial condition Fit II given in Table 2) over 10 5 yr. It appears 
to be rigorously stable which according to the theory of the 
fast indicator is indicated by MEGNO rapidly converging to 
the value of 2. Because the low-mass innermost planet re- 
volves very close to the parent star, its influence on the secu- 
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lar dynamics is negligible. Thus in the next experiments we 
considered the three- Jovian system only. The results of the 
integrations conducted over 50 Myr are shown in Fig. [7] Such 
period of time corresponds to ~ 5 • 10 6 P C and is long enough 
to detect possible secular instability. Still the MEGNO con- 
verges perfectly to 2 and that indicates a stable quasi-periodic 
solution. Indeed, the eccentricities and semi-major axes are 
almost constant over this time. 

In order to illustrate the dynamical environment of the 

Ara system, we also computed the dynamical maps in 
terms of the Spectral Number and maxe in the (at,, em- 
plane (Fig. [8}. These maps reveal that the nominal posi- 
tion of planet b is close to the island of 2:1 MMR with the 
planet e. Yet none of the critical angles of this resonance li- 
brate in the nominal best-fit configuration. We found only 
some signs of the apsidal anti-corotation (see the left-lower 
panel in Fig.0. Such dynamical character of the /j Ara sys- 
tem is puzzling in the light of what is known about four ex- 
trasolar systems presumably locke d in the exact 2:1 MMR, 
i.e., Gliese 876 (iRivera et alJl2005l) . H P 82943 dMavor et all 
2004; Lee et alJ 120061) HP 128311 dVo gt et alJ 120051) and 
HP 73526 (Tinnev et al. 2006). The /j Ara system seems only 
close to the border of this resonance. A possible explanation 
of this behavior is given in the next section. 

One should be aware that K c (and the resulting minimal 
masses) as well as P c can be significantly varied within the 
la confidence range of the best fit. Moreover, not all Kep- 
lerian fits shown in Fig. |5] are necessarily stable. Thus the 
GAMP search would be very helpful to get a more detailed 
self-consistent statistics of the stable solutions. In this work 
we skipped such a test due to its significant numerical expense 
caused by the extremely different orbital periods. Instead, we 
carried out a more direct check of the stability preserving the 
full accuracy of the four-planet fits. First, the ~ 1000 different 
best Keplerian fits gathered in the hybrid search were refined 
by the Newtonian model of the RV. Next, we performed long- 
term MEGNO integrations of these self-consistent solutions. 
At this stage, the planet d has been skipped. The integration 
time t m ~ 3 • 10 5 f c (<~ 3 Myr) is long enough to detect un- 
stable solutions related to the short-term two- and three-body 
MMRs and also strongest secular resonances. A detection of 
all secular instabilities would require much longer integration 
times ~ 10 7 — 10 8 yr because the secular periods, as the apsidal 
period of the outermost orbit, are ~ 10 5 yr. 

The results of this experiment are shown in the top row of 
Fig- ED The best-fit solutions have an rms ~ 2.27 m/s. Glob- 
ally there is only a little improvement of the Newtonian fits 
when compared to the Keplerian solutions. The mutual inter- 
actions are not yet clearly evident in the AAT data. The over- 
all distribution of the Newtonian fits in the (a c , e c ) -plane mim- 
ics the projection of the kinematic solutions onto the (P c , em- 
plane (see Fig. There is a well defined minimum of rms 

[and (Xv) 1 ^ 2 ] m me (fle,bie e ,b)-pl anes , but the semi-major axis 
of planet c may be varied over 2 AU within the limit of rms 
< 2.3 m/s. We examined the stability of all such fits. The sta- 
ble initial conditions (with \ (Y)(t m ) — 2| < 0.001) are marked 
with filled yellow circles. We notice that the stable fits appear 
for e e ,b < 0.1. It remains likely that by altering the orbital 
phases or other parameters (like masses) within the accept- 
able error bounds, we can find reasonably precise and stable 
solutions also for <? e ,b > 0.1. However, this procedure would 
require the self-consistent GAMP-like search. 



5. FOUR-PLANET MODEL REVISITED 

Shortly after submitting the original manuscript, we learned 
about an independent work by Pepe et all (120061) . Our conclu- 
sions are in a great agreement with the results of their analy- 
sis although these authors study a completely different and 
independe nt set of the RV data: the measurements published 
bv lMcCarthv et all (120041) . data gathered with CORALIE and 
observations collected during ~ 2 year campaign with the 
HARPS spectrometer. Having access to the new extended ob- 
servations, we can discuss some of the conclusions derived on 
the basis of the new AAT data alone. 

In a number of experiments , we co nsidered three data sets: 
SI — 108 points by the AAT ( iButler et alJl2006l) as analyzed 
already, S2 — the set SI extended by the CORALIE mea- 
surements (with errors rescaled by Oj ^3.5 m/s) and 78 mea- 
surements from HARPS (the original errors are unmodified) 
and S3 — the set composed of only the AAT and HARPS 
data. Note, that from the CORALIE and HARPS observa- 
tions (Pepe et all 120061) . we substracted the mean of all RV 
measurements in the given set, respectively. 

Using the most precise S3 set which covers the whole ob- 
servational window, we carried out the hybrid search for the 
four-planet Keplerian fits including two independent instru- 
mental RV offsets. As in the previous test, the orbital periods 
of all planets are bounded to the range of [8,6500] days. 
The best Keplerian fit is similar to Fit II, i.e., in terms of the 
parameter tuples of Eq. 1, (3.118 m/s, 9.636 days, 0.164, 
21 1.68 deg, 302.852 days), (10.849 m/s, 307.937 days, 0.127, 
180.752 deg, 3323.904 days), (37.147 m/s, 649.500 days, 
0.000, 231.352 deg, 2734.35 days), (24.056 m/s, 
4293.506 days, 0.054, 90.829 deg, 3319.1 days), for planets 
d,e,b,c, respectively and velocity offsets Vq = —10.932 m/s, 
Vi = 1 1 .209 m/s (r p are shifted by JP2,450,000). This fit has 

(Xv) 1 ^ 2 ~ I- 48 and an rms ~ 2.51 m/s. 

Next, we refined ~ 2500 different solutions gathered in the 
hybrid search using the Newtonian code. Similarly to the case 
of the AAT data only (SI), we also performed the stability 
check of the best-fit configurations. For the S2 set (see the 
middle panels in Fig. [9}, the smallest rms of the best-fit solu- 
tions is ~ 3.66 m/s, significantly more than that one obtained 
for the S 1 set alone (2.3 m/s); likely due to a relatively low ac- 
curacy and scatter of the CORLIE data. The best fits to the S3 
data yield an rms ~ 2.5-2.7 m/s also larger by ~ 0.3-0.4 m/s 
from the rms obtained for the S 1 set. Also their (x 2 ) 1//2 ~ 1 -5 
suggest that the errors of the HARPS measurements could 
be underestimated in the long-run. Hence finally we mini- 
mized the rms instead of (x 2 ) 1//2 - In this case (see the bottom 
panels of Fig. |9j, the rms of the best fits is smaller than for 
the SI set alone, ~ 2.2 m/s. For instance, a stable solution, 
given in terms of tuples (m [Mj],a [AU],e,(0 [deg],5Vf [deg]): 
(0.032, 0.09286, 0.135, 204.663, 272.050), (0.506, 0.9393, 
0.083, 204.913, 293.54), (1.696, 1.533, 0.062, 45.56, 4.039), 
(1.926, 5.134, 0.053, 103.06, 156.057), for planets d, e, b, 
c, respectively, with velocity offsets Vq = —8.466 m/s and 
V\ = 13.695 m/s yields an rms ~ 2.22 m/s. 

Still, for all the analyzed sets SI, S2 and S3, the orbital pa- 
rameters (a c ,e c ) are allowed to vary in relatively wide ranges 
qualitatively similar to the ones we determined for the AAT 
data only. Yet the shapes of the minima are much more ir- 
regular than those found for the set S 1 . Clear limits of stable 
solutions in the (ab,e,eb,e)-planes are evident. 

Having the ensemble of stable fits, we also try to explain the 
apparent lack of locking of the subsystem e-b in the 2: 1 MMR. 
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Because the elements of the inner giants are well constrained, 
we selected three best fits to the S3 data (all yielding an rms 
~ 2.2-2.3 m/s) with a c ~ 4.7,5,6.3 AU, respectively; also 
Fit II (a c ^5.5 AU) may be put in this short sequence. Next, 
we calculated the dynamical maps in the (a^etO-plane m the 
neighborhood of the selected solutions. These maps (Fig. 1101 
reveal very sharp borders of stable motions around et,,e ~ 0.1, 
in an excellent agreement with the MEGNO tests (Fig.|9}- The 
outermost planet strongly modifies the dynamics of the e-b 
pair. For a moderate distance of planet c from the e-b subsys- 
tem, there is no evident 2:1 MMR island at all. It shows up for 
a c larger than 5-5.5 AU and for relatively large eccentricities 
(see also Fig.|5J- However, such large eccentricities are ruled 
out in the model by the stability constraints. 

5.1. The le:lb MMR hypothesis 

Remarkably, the Keplerian and A^-body fits to all the data 
sets (S1,S2,S3) reveal a number of low rms solutions corre- 
sponding to the lb:le MMR of the planet s e and b. This pos- 
sibility is also investigated by Peo e et all (120 06) but they did 
not report any stable 1:1 MMR solutions. 

The 1 : 1 MMR configurations are quite frequent in the So- 
lar system. There are also speculati ons on the existence 
of such extrasolar configurations tL aughlin & Chambers 
1 200 11 iNau enberg M] 120021 Idozdziewski & Konackil 120061 
Ford & Gaudi 2006). For a closer analysis we selected the fits 
to the S3 set. The Newtonian best-fit solutions with a e ~ a\, 
are illustrated in the top-left panel of Fig.[^ Their quality is 
comparable with the ones related to the 2e:lb MMR. Never- 
theless, all these le:lb MMR configurations are strongly un- 
stable leading to a quick collision between the planets. Due 
to extremely strong dynamical interactions, even small errors 
of the phases or other parameters may be critical for the sys- 
tem stability. Thus in a close proximity of apparently unsta- 
ble solutions, stable systems consistent with the RV observa- 
tions may still exist. We tried to "stabilize" such fits with 
GAMR in terms of the three- Jovian planets model (that is 
again due to the efficiency reasons). Unfortunately, we did 
not find any long-term stable best-fits as precise as the ones 
of the 2e:lb MMR configurations. However, there exist sta- 
ble le:lb MMR solutions with reasonably small rms ~ 5 m/s 
(~ 4.5 m/s for the whole, four-planet system). A synthetic 
curve of an example configuration fitting well the most pre- 
cise AAT and HARPS measurements, is illustrated in the top- 
right panel of Fig. (the orbital elements are given in the 
caption). A striking result is illustrated in the dynamical map 
of that best-fit (see the right-bottom panel of Fig. lilt . It re- 
veals the island of stable motions extending over 0.2 AU in 
a e and covering the entire range of e e \ The selected solution 
is stable over the secular time-scale. We show the results of 
1 Gyr direct integration of the system including the planets 
b,c and e (see the left-bottom panel in Fig. II II . The apsides 
of e-b subsystem are anti-aligned, librating around 180° with 
the semi-amplitude of ~ 70° and with et,,e varying up to 0.5 
and <? c ~ 0.1. We also found other solutions with the apsides 
librating about centers different from 180°. The le:lb MMR 
island has a very complex dynamical structure, particularly 
in the vicinity of the best fit. Direct integrations show that 
weakly chaotic solutions in that zone may lead to a sudden 
collision of planets after hundreds Myr of apparently stable 
orbital behavior. Let us note that to refine the stable fit, the 
MEGNO was integrated over ~ 10 5 P c . 

Still, there is an open question whether the fit quality should 



be the final argument to rule out the lb:le MMR configura- 
tions. We think that the problem deserves a further study. Let 
us recall that we searched only for coplanar solutions. The 
stabilit y may be easily preserve d for inclined configurations 
(Gozdziewski & Konacki 2006). Yet the system architecture 
involving planets b and e in the 1 : 1 MMR may still permit the 
existence of other bodies besides the well established planet c. 

6. CONCLUSIONS 

The long-term radial velocity surveys of the Sun-like 
stars constantly reveal more and more exciting features 
of the planetary systems. The /j Ara system may be 
the second known fo ur-planet configuration, after 55 Cnc 
(McArthuret al. 2004). Remarkably, in such a multi-planet 
system the orbits are close to circul ar ones, sim ilarly to three- 
planet systems of H P 37124 (IVogt et all2005l) and HD 68930 
tLovis etafll2006l) resembling the Solar system architecture. 
The alternative best-fit three-planet configurations may con- 
tain two Jupiter-like planets in the 4: 1 MMR. In that case, the 
eccentric orbits (eb,c ~ 0.3) would be localized in a dynami- 
cally active region of the phase space; in fact, on the edge of 
dynamical stability. Besides the worse fit quality, it could be 
a heuristic argument against the three-planet system. Obvi- 
ously, the key for the proper understanding of the orbital ar- 
chitecture is the improved precision of observations and their 
extended time span. Curiously, the new data reveal a Jovian 
planet that has the orbital period two times shorter than the 
companion already detected a few years ago. It remains an 
open question whether the two inner companions (the plan- 
ets b and e) in the four-planet configuration are locked in the 
2:1 MMR. Most likely, the presence of the massive compan- 
ion c prevents the creation of the 2:1 MMRs island known 
in the other four extrasolar systems (Gliese 876, HD 82943, 
HD 128311, and HD 73526) presumably locked in this res- 
onance. In any case, even the proximity of their orbits to 
such particular dynamical state can be counted as the fifth 
occurrence of the 2:1 MMR among ~ 20 multi -planet sys- 
tems known up to date. It may indicate a universal dynam- 
ical mechanism governing the creation and orbital evolution 
of extrasolar planetary systems. We also found stable config- 
urations related to the 1 : 1 MMR of the inner Jovian planets 
in eccentric orbits. However, such fits are significantly worse 
than those derived for the system with quasi-circular orbits. 
The results of our analysis permit us to conclude that a few 
years of observations are still required to constrain the outer- 
most orbit without a doubt. Lookin g at the orbits o f the fj Ara 
planets and recalling the results of iLaskan (12000 ). we see a 
free space in the range of ~ (0.2,0.8) AU in which new plan- 
etary objects may yet exist. Extensive numerical simulations 
concerning the less constrained parameters of the outermost 
planet would be necessary to answer this question. 
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TABLE 1 

Primary parameters (K,P,e,(o,T P ) of the three-plane t Keplerian model (i n Eq. 1) and the inferred astmcentric orbital elements of 

THE BEST-FIT I FOUND IN THIS PAPER FOR THE RV OF (1 ARA jBUTLER ET AL.I2006I) . THE JITTER ESTIMATE CJ = 3.5 M/S, THE MEAN MEASUREMENT 
ERRORS (cj) = 1.87 M/S, AND THE MEAN COMBINED DATA ERROR (OJ 2 + {a} 2 ) 1 / 2 = 4.02 M/S. NUMBERS IN PARENTHESES ARE FOR THE lo" ERRORS 
ESTIMATED ON THE BASIS OF THE BEST FIT STATISTICS GATHERED IN THE HYBRID SEARCH (SEE THE TEXT AND FlG.0. BY DEFINITION, T e AND CO 
ARE NOT WELL CONSTRAINED WHEN e IS SMALL. INSTEAD, THE ERROR OF THE ORBITAL PHASE, OR THE MEAN LONGITUDE, X(t ) = W (/ ) + CO, 
WHERE <M (f ) IS THE MEAN ANOMALY, IS GIVEN. THE EPOCH OF THE FIRST OBSERVATIO N fp = JD 2,45 1 , 1 1 8.8 9. THE REFERENCE EPOCH T IS 

JD 2,440,000. Mass of the parent starM* = 1.15 M @ JButler et al.I2006I) . 



Best Fit I HD160691b HD160691c HD160691d 



P [days] 


632.013 (6) 


2544.47 (60) 


9.6369 (0.005) 


K [mis] 


37.97(1) 


17.67 (1) 


3.03 (0.50) 


e 


0.260 (0.07) 


0.471 (0.05) 


0.236 (0.2) 


CO [deg] 


259.46 


183.63 


304.20 


T p [JD-7/ ] 


12137.86 


13536.77 


10545.71 


X(to) [deg] 


39.05 (7) 


201.54 (7) 


115.54 (9) 


msini [mj] 


1.70 


1.15 


0.034 


a [AU] 


1.51 


3.80 


0.09285 


ixl) 1 ' 2 




1.01 




rms [m/s] 




3.98 




Vq [m/s] 




-0.638 (0.8) 
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TABLE 2 

Primary parameters (K,P,e,a,T P ) of the four-planet Keplerian best fit II (Eq. 1) and the inferred astrocentric orbital elements. 
This fit is dynamically stable. Numbers in parentheses are for the la errors estimated on the basis of the best fit statistics 
gathered in the hybrid search (Fig.|5J. See caption to Table 1, Fig.[5]and the text for more details. 



Best Fit II 


HD 16069 lb 


HD160691c 


HD 16069 Id 


HD160691e 


P [days] 


646.485 (1.5) 


4472.967 (1300) 1 


J.6369 (0.005) 


307.475 (1.5) 


/C" [m/s] 


35.871 (1) 


27.178 (10) 


2.826(1) 


13.195 (2) 


e 


0.0001 (0.05) 


0.027 (0.12) 


0.184 (0.2) 


0.079 (0.06) 


co [deg] 


223.003 


154.065 


314.050 


252.624 


T v [JD-r ] 


12721.839 


14171.256 


10632.575 


13070.393 


X(f ) [deg] 


50.386 (2) 


268.400 (36) 


121.225 (10) 


127.752(10) 


msirn [mj] 


1.677 


2.423 


0.032 


0.480 


a [AU] 


1.535 


5.543 


0.09286 


0.934 


(X?) I/2 




0.627 






rms [m/s] 




2.276 






Vo [m/s] 




-13.069 (10) 





Fig. 3. — The dynamical maps in the (a c ,e c )-plane for the three-planet Keplerian model of the /j Ara system. The large crossed 
circle marks the parameters of the best Fit I. The left panel is for the Spectral Number, \ogSN. Colors used in the logSN map 
classify the orbits — black indicates quasi-periodic regular configurations while yellow strongly chaotic ones. The right panel 
marked with maxe e is for the maximal eccentricity of planet e attained during the integration of the system. The thin line marks 
the collision curve for planets b and c, as determined by a\,(l +eb) = fl c(l — e c)- The low-order MMRs of planets b and c are 
labeled. The integrations are conducted over ~ 10 4 P C . The resolution is 400 x 120 data points. 

Fig. 4. — The dynamical maps in the (a c , e c )-plane and the parameters of the best-fits derived by GAMP. The edge-on, two-planet 
model is assumed. The nominal initial condition is marked by a large crossed circle. Its osculating elements {m p ,a,e,(0, M ) at 
the epoch of the first observation are: (1.572 m h 1.514 AU, 0.251, 253.222 deg, 145.937 deg) and (1.182 mj, 3.858 AU, 0.332, 
189.385 deg, 17.771 deg), for the inner and outer planet, respectively, and Vo = —0.85 m/s. This solution yields (xl) 1 ^ 2 ~ 1.17 
and an rms ~ 4.7 m/s. The relevant parameters of other fits within the formal la confidence interval [(x 2 ) 1//2 < 1.18 and an rms 
4.8 m/s] of the initial condition are marked by small circles. See the caption to Fig.[3]for an additional explanation. 

Fig. 5. — The parameters of the best-fit solutions to the four-planet Keplerian model of the RV of /j Ara, projected onto the (P,e)- 
and (P,/T)-plane. The values of (x 2 ) 1//2 °f tne best-fit solutions are marked by the size of symbols (smaller {% 2 ) l ^ 2 2 — larger 
circles). The largest circle is for (x 2 ) 1//2 equal to 0.626; smaller symbols are for la solutions with (x 2 ) 1//2 < 0.645, 2a solutions 
with (Xv) 1//2 < 0.68, and 3a solutions with (%l) 1 ^ 2 < 0.73 (smallest circles), respectively. The elements of the best fit found in 
the search are marked by two crossing lines. 

Fig. 6. — The synthetic signal of the four-planet Keplerian best Fit II, see Table 2. Subsequent panels (from the top) are for the 
synthetic RV signal of the four-planet model, the period-phased RV signals of the innermost, Neptune-like companion d, the new 
companio n e, the most m assive planet b, and the outermost planet c, respectively. The open circles are for the RV measurements 
from|Butler et al. (2006). The error bars include the internal errors added in quadrature to stellar jitter of 3.5 m/s. The next panel 
is for the Lomb-Scargle periodogram of the RV data in the range of the short periods. The last two panels are for the periodogram 
of the residual signal after subtracting the contribution of Jovian planets, and the histogram of (% 2 ) 1 / 2 derived in the test of 
scrambled residuals. 

Fig. 7. — Evolution of MEGNO (the top-left panel labeled by (Y)), and orbital elements of the three-planet configuration 
described by Fit II given in Table 2 (only Jovian planets are considered). A perfect convergence of MEGNO up to ~ 50 Myr 
indicates a rigorously stable solution over secular time scale. The subsequent panels are for the eccentricities, the angle measuring 
the apsidal anti-alignment of orbits e and b, and the semi-major axes. 

Fig. 8. — The dynamical maps in the (at,,<?b)-plane of the /j Ara system for the four-planet best Fit II (Table 2). The thin line 
marks the collision curve for planets b and e. See the caption to Fig.[5]for an additional explanation of the plots. 

Fig. 9. — The statistics of the best, self-consistent Newtonian fits and their dynamical stability. The subsequent panels are for 
the projections of the best fit, osculating elements of the Jovian planets at the epoch of the relevant first observation. The quality 
of fits, expressed by their rms, is marked by the size of symbols (circles) and labeled in the plots. White (yellow-in the color 
version of the figure) circles are for dynamically stable best-fit solutions. The orbital stability is examined through MEGNO 
integrations over ~ 3 • 10 5 P C , for all solutions with lowest rms (marked in the panels by largest circles and labeled accordingly). 
Best fit-solutions marked as stable have | (Y) — 2| < 0.001 at the end of the integration period. The upper row is for the S 1 data set 
[the AAT measurements published by (Butl er et alJ l2006)l. The middle row is for the S2 data set (AAT+CORALIE+H ARPS ) . 
The bottom row is for the S3 data set (AAT+HARPS); note, that in that case the rms was minimized instead of (x 2 ,) 1 / 2 - Crossed 
lines marks the elements of the fits with lowest rms. 
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Fig. 10. — The dynamical maps in the (a e ,e e ) -plane computed for the best fits to the S3 data set (including AAT and HARPS 
observations, see Fig. 9). The maps are computed for the following osculating elements of the Jovian planets, given in terms of 
tuples (m [Mi], a [AU],e,co [deg],M [deg]): the left panel is for (0.430, 0.937, 0.086, 205.973, 282.203), (1.692, 1.532, 0.015, 
314.000, 95.498), (1.704, 4.702, 0.088, 161.445, 81.240); the middle panel is for (0.549, 0.940, 0.077, 206.297, 298.451), (1.709, 
1.534, 0.101, 47.902, 1.628), (1.808, 4.969, 0.073, 128.199, 126.093); and the right panel is for (0.485, 0.939, 0.085, 205.241, 
289.819), (1.686, 1.533, 0.042, 41.633, 8.237), (2.643, 6.306, 0.122, 23.664, 263.507), respectively. The thin line marks the 
collision line for planets b and e. See the caption to Fig.[3]for an additional explanation of the plots. 

Fig. 1 1 . — The dynamical analysis of the orbital configuration of the /j Ara system involving planets e-b in 1 : 1 MMR. The top-left 
panel is for the projections of the best fits at the (<?b, a e )-plane of osculating elements, as derived for the S3 data set (see the text for 
an explanation). The top-right panel is for the synthetic curve of a stable solution, refined by GAMP integrations over ~ 10 5 P C . 
The osculating elements of the three-Jovian system are: (1.196, 1.52896, 0.3195, 180.208, 213.185), (0.602, 1.49056, 0.305, 
31 1.601, 69.827), (1.939, 5.035, 0.016, 136.75, 110.0), respectively, given in terms of tuples (m [Mi], a [AU],e,co [deg],5Vf [deg]), 
at the epoch of the first observation, and the offsets are V\ — —11.395 m/s, V% — 11.124 m/s. The bottom-left panel is for the 
time-evolution of 9 = G3 e — EtSb, in the best-fit configuration, during 1 Gyr. The bottom-right panel is for the dynamical map in the 
vicinity of the best fit; its position is marked by the crossed circle. 
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